Exploring genetic and expression differences between physiologically extreme ecotypes: comparative genomic hybridization and gene expression studies of Kas1 and Tsu1 accessions of Arabidopsis thalianaT. E. Juenger etal.Genomic studies of Arabidopsis accessions

Recent studies have documented remarkable genetic variation among Arabidopsis thaliana accessions collected from diverse habitats. Of particular interest are accessions with putatively locally adapted phenotypes – that is, accessions with attributes that are likely adaptive at their sites of origin. These genotypes may provide insight into the genetic basis of adaptive evolution as well as allow the discovery of genes of ecological importance. We studied the physiology, genome content and gene expression of two physiologically extreme accessions (Tsu-1 from Tsushima, Japan and Kas-1 from Kashmir, India). Our study was conducted under two levels of soil moisture and accompanied by physiological measurements to characterize early responses to soil drying. Genomic hybridizations identified 42 503 single feature polymorphisms (SFP) between accessions, providing an initial screen for genetic differences. Transcript profiling identified a large number (5996) of genes exhibiting constitutive differences in expression including genes involved in many biological pathways. Mild soil drying resulted in only subtle physiological responses but resulted in gene expression changes in hundreds of transcripts, including 352 genes exhibiting differential responses between accessions. Our results highlight the value of genomic studies of natural accessions as well as identify a number of candidate genes underlying physiological differences between Tsu-1 and Kas-1. Key-words: Acclimation; affymetrix; drought; expression polymorphism; hybridization polymorphism; soil moisture; water-use efficiency. INTRODUCTION As sessile organisms, plants are constantly challenged by a variety of abiotic stresses from which they cannot flee. As such, they are a fascinating system for studying the physiological mechanisms underlying abiotic stress tolerance, including both acclimation responses and local adaptation. An especially critical environmental stress for plants is the availability of soil moisture, as water availability is fundamental to almost all aspects of plant physiology. Plant water status has a strong and direct impact on C3 photosynthesis and carbon fixation through stomatal regulation of gas exchange, as well as important consequences for plant growth, phenology and susceptibility to other abiotic and biotic stresses (Rizhsky et al. 2003). More broadly, water availability and atmospheric demand, interacting with temperature, are fundamental determinants of plant distribution abundance, and productivity worldwide (Walter 1964; Whittaker 1975). As a result, crop yields are commonly reduced by water limitations to less than half of potential yields (Boyer 1982; Gleick 1998). Therefore, research exploring the molecular genetic and physiological basis of tolerance to soil drying is critical to understanding and improving plant function under stressful conditions. The model angiosperm Arabidopsis thaliana has emerged as a valuable tool in deciphering how plants respond to environmental stresses, including drought and desiccation (Zhang, Creelman & Zhu 2004). Decades of research has shown that cell signalling and gene expression networks underlie physiological adjustment and are likely to affect plant performance under stress. As an example, a host of studies in Arabidopsis have documented gene expression changes induced by dehydration-related treatments. Two cis-acting DNA sequence elements clearly contribute to water deficit-induced gene expression in the Columbia accession, suggesting that regulatory evolution may be a key process in drought adaptation in Arabidopsis. The Correspondence: T. E. Juenger. E-mail: [email protected] Plant, Cell and Environment (2010) doi: 10.1111/j.1365-3040.2010.02146.x © 2010 Blackwell Publishing Ltd 1 abscisic acid (ABA)-response element (ABRE) is important for ABA-dependent changes in gene expression (Uno et al. 2000; Bray 2004), and many ABA-responsive genes are likely to play a role in drought acclimation. Moreover,ABA functions in stomatal regulation, control of growth and osmolyte accumulation (Verslues & Zhu 2007). The dehydration-response element/C-repeat (DRE/CRT) is essential for ABA-independent induction of many desiccation-responsive genes (Yamaguchi-Shinozaki & Shinozaki 1994). Several of the transcription factors (DREB/ CBF) involved in DRE/CRT-responsive gene expression have been cloned and a number of target genes characterized (Shinozaki & Yamaguchi-Shinozaki 2000; Seki et al. 2001, 2002; Bray 2002a; Shinozaki et al. 2003). Genetic manipulation of drought-responsive transcription factors and/or their downstream targets has resulted in increased plant performance under specific stress treatment conditions (reviewed in Umezawa et al. 2006). Despite this progress, we still have much to learn. Although the existing literature on stress induction has been critical for benchmarking gene function (Bray 2002a,b, 2004; Kilian et al. 2007), lab treatments (desiccation of cut leaves on benchtop environments), short response times (1–24 h) and limited genetic diversity (wild-type, knockout or overexpression lines) may limit the scope of inference from these experiments. In particular, it is possible that many of the effects observed are better considered ‘stress shock’, given the timescales of stress imposition and the lack of acclimation. It is likely that ‘shock’ studies have touched only the surface of how plants respond to diverse abiotic stresses.As such, we have initiated a number of studies implementing more realistic progressive soildrying manipulations and utilizing a diversity of accessions. A. thaliana provides a unique opportunity to explore adaptive evolutionary responses to drought as it has an extensive geographical distribution and has experienced a wide-range of climatic selective regimes for thousands of generations (McKay, Richards & Mitchell-Olds 2003; Koornneef, Alonso-Blanco & Vreugdenhil 2004; Bouchabke et al. 2008). We have completed several common garden and quantitative genetic and quantitative trait loci (QTL) mapping experiments focused on plant–water relations and integrative water-use efficiency measures in Arabidopsis (McKay et al. 2003; Hausmann et al. 2005; Juenger et al. 2005, 2006; Christman et al. 2008; McKay et al. 2008). These studies have characterized the general range of variability as well as identified a number of accessions with extreme physiological traits. We have centred our recent efforts on the Kas-1 and Tsu-1 accessions, as they represent the highest and lowest identified water-use efficiencies in diversity panels, respectively (McKay et al. 2003), and the climate of the accession sites of origin differs greatly in both precipitation and temperature. The site of origin of Tsu (Tsushima, Japan) has high water availability throughout the growing season and the Kas site of origin (Kashmir, India) has very limited precipitation inputs during the growing season (McKay et al. 2008). To provide a tool to understand the genetic basis of this putative adaptive differentiation, Kas-1 and Tsu-1 were reciprocally crossed to create a new recombinant-inbred mapping population (McKay et al. 2008).The Kas-1 ¥ Tsu-1 mapping population provides a powerful new resource for mapping QTL underlying natural variation and for dissecting the genetic basis of water-use efficiency differences. Here, we extend previous studies of drought-induced gene expression to incorporate more realistic stress treatments and putatively locally adapted plant material. We utilize genomic hybridization and gene expression approaches using Affymetrix GeneChip microarrays (Santa Clara, CA, USA) and the Tsu-1 and Kas-1 accessions. Our goals are to characterize natural physiological responses to soil drying in Arabidopsis, identify putative sequence differences between physiologically extreme accessions and explore patterns of transcript responses to initial and mild stress conditions. MATERIALS AND METHODS Plant growth and drydown Tsu-1 (CS 1640) and Kas-1 (CS 903) were grown in 6 ¥ 6 ¥ 5 cm plastic pots filled with ~150 mL Profile porous ceramic rooting media (Profile Products LLC, Buffalo Grove, IL, USA). Prior to planting, dry weight of each pot and soil (DWpot & soil) was measured and the media was saturated with a complete nutrient solution (Epstein & Bloom 2005) for several days. Saturated pots were covered and allowed to drain by gravity for >24 h until a constant field capacity weight (FCpot & soil) was obtained for each pot. Three seeds were planted in each pot and dark stratified at 4 °C for 3 days. Germination was on 9–10 August 2006. Plants were thinned to one per pot and were grown and treatments were applied in a controlled environment chamber with 10 h photoperiod of 330 mmol m-2 s-1 photosynthetic photon flux density. Daytime/night-time temperatures were 23/18 °C with 60% relative humidity. During growth, all plants were fertilized with the nutrient solution every two days and watered daily.The drydown began on 26 August 2006 with all plants in rosette stage. Plants were grown in three replicate blocks with treatments and accessions randomly located within each block. There were six subsample plants each of Tsu and Kas in each treatment in each block (72 plants total). During the drydown treatment, each pot was weighed daily at ‘pre-dawn’ (WTpot & soil) and distilled water was added by pipette at the base of the rosette to bring the percent water remaining (WatRem%) in each pot to the target level for that day (e.g. 90% on day 1, 80% on day 2, etc.) until a final drying treatment level of ~50% was reached by all pots on the fifth day.A major goal of the slow drydown treatment was to allow the intact, growing plants to acclimate to the imposed soil moisture stress. WatRem% was calculated as 100 ¥ (WTpot & soil – DWpot & soil)/ (FCpot & soil – DWpot & soil). Wet (control) pots were treated identically, except WatRem% was maintained near 90% for the entire drydown period. At ‘dusk’, prior to the 2 T. E. Juenger et al. © 2010 Blackwell Publishing Ltd, Plant, Cell and Environment harvest day, plants were enclosed in covered trays to minimize evapotranspiration. Covered trays were transported to the lab the following morning, final pot and soil weights were taken, WatRem% at harvest were calculated and plants were harvested. Harvest occurred during the first several hours of the normal daily light period. Harvest and physiological measurements Plants of both lines and both treatments in each block were harvested together. Rosettes were excised and fresh weight of rosette and any bolt stem were determined separately. Immediately, ~80 mg (~40% of the leaves on the rosette) of healthy, fully expanded leaves were detached from the rosette and placed in RNAlater (Ambion, Inc., Austin, TX, USA), several additional leaves were removed and placed in a psychrometer chamber for water potential measurement (below) and 2–3 leaves were removed for water content (WC) and relative water content (RWC) measurements. For WC and RWC measurements, leaf fresh weight was immediately determined on a microbalance, leaf bases were placed in distilled water in a microfuge vial and the vial and leaves were enclosed in a larger sealed vial to allow leaves to reach full hydration in the dark for ~18 h. Turgid weight of the hydrated leaves was then determined, leaves were dried and dry weight determined using the same microbalance. WC was calculated as 100 ¥ (fresh weight dry weight)/(dry weight) and RWC as 100 ¥ (fresh weight dry weight)/(turgid weight dry weight) (Boyer 1995). Remaining leaves were weighed fresh and after drying so that total rosette and total leaf weight, both fresh and dry, could be calculated by summing values from material used for each type of measurement. After drying, leaves were used for C and N content and stable isotope analyses at the UC Davis Stable Isotope Facility (http:// stableisotopefacility.ucdavis.edu). Carbon isotope composition (dC) is given relative to the PeeDee Belemnite standard and composition is used rather than discrimination (D) because the isotopic composition of carbon dioxide in the ambient air, which is required to compute D, was extremely variable (see McKay et al. 2003). Leaf water potential (Ytot) was measured using the excised leaves, individually calibrated (Brown & Bartos 1982) thermocouple psychrometers (Merrill Specialty Equip., Logan, UT, USA) and stainless steel chambers (see Donovan, Linton & Richards 2001; Donovan, Richards & Linton 2003). Entire leaves were placed in chambers within 30 s of excision and chambers were sealed and suspended in a water bath to minimize temperature gradients during measurement. Mature leaf tissue was used to minimize any growth effects on leaf Ytot (Boyer 1995). Psychrometer outputs were logged hourly (CR7 data logger; Campbell Scientific, Logan, UT, USA) and leaf Ytot was determined after equilibration (~24 h). After Ytot was measured, the sealed chambers were immersed in liquid nitrogen to rupture cell membranes of the leaves and an estimate of bulk leaf osmotic potential (Ysol) was determined during a second 24 h equilibration period. For each subsample plant, leaf turgor (Yp) was calculated as Ytot Ysol. Similar to leaf measurements, soil water potential (Ysoil) was determined on a sample of soil, with roots, taken from the centre of a subset of the pots and placed in stainless steel chambers with thermocouple psychrometers. All water potentials were calculated from psychrometer mV outputs using the Brown & Bartos (1982) model that accounts for zero offset (always between 0.2 and 0.2 mV) and temperature. Genomic DNA hybridizations to ATH1 array Comparative genomic hybridizations were completed following the protocols outlined in Borevitz (2006) and using the Invitrogen (Carslbad, CA, USA) BioPrime Labeling System. Genomic DNA was extracted from seedlings of Tsu-1, Kas-1 and Col-0 (CS 60000) accessions using the Qiagen DNAeasy (Germantown, MD, USA) plant kit. Approximately 300 ng of genomic DNA was added, on ice, to 60 mL of 2.5 ¥ random primer solution and made to a final volume of 132 mL with distilled water. The mixture was denatured by incubation at 99 °C for 10 min and immediately placed on ice for 5 min. Next, 15 mL of 10 ¥ deoxyribonucleotide triphosphate solution (with biotin-deoxycytidine triphosphate) and 3 mL Klenow were added to the denatured DNA mixture. The reaction was incubated in a PCR block at 25 °C for 16 h and the reaction was terminated by the addition of 15 mL stop solution. Labelled DNA was precipitated by the addition of 20 mL of 3 mol L-1 sodium acetate and 400 mL cold 100% EtOH. This solution was incubated on ice for 2 h and centrifuged at 15 000 g for 10 min followed by a wash with 70% EtOH. The pelleted DNA was dried and resuspended in deionized distilled H20. The labelled DNA was then treated the same as labelled cRNA using standard hybridization protocols for gene expression (see methods below). Genomic hybridizations to the Affymetrix ATH1 array were completed for six samples of each accession for a total of 18 arrays. Hybridizations were completed in two separate batches of nine arrays. The genomic hybridization data from this study are accessible through GEO Series accession number GSE20340 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE20340). Gene expression studies We used a subset of plants from the drydown experiment in our studies of soil water deficit-regulated gene expression. In brief, our experimental design was a fully factorial experiment involving two accession (Tsu-1 and Kas-1), two treatments (control and soil drying), with six biological replicates (2 ¥ 2 ¥ 6 = 24 arrays) from a single experimental block. As described above, fully expanded rosette leaves were sampled on RNAlater (Ambion, Inc.) and total RNA was extracted using Qiagen RNAeasy kits. Samples for mRNA profiling studies were processed by Asuragen, Inc. (Austin, TX, USA) according to the company’s standard Genomic studies of Arabidopsis accessions 3 © 2010 Blackwell Publishing Ltd, Plant, Cell and Environment operating procedures.The integrity of total RNA was qualified by Agilent Bioanalyzer 2100 capillary electrophoresis (Palo Alto Santa Clara, CA, USA) and used for preparation of biotin-labelled targets (cRNA) using a MessageAmp II-based protocol (Ambion Inc.). The cRNA yields were quantified by ultraviolet spectrophotometry and the distribution of transcript sizes was assessed using the Agilent Bioanalyzer 2100 capillary electrophoresis system. Labelled cRNA was fragmented and used for array hybridization and washing, according to the standard Affymetrix protocol. In brief, labelled cRNA was resuspended in 5¥ fragmentation buffer and incubated at 94 °C for 35 min then stored on ice. The hybridization cocktail and the fragmented cRNA mixture were heated to 99 °C for 5 min, and incubated at 45 °C for 5 min. After a final spin to collect the samples, hybridization to arrays was carried out at 45 °C for 16 h in an Affymetrix Model 640 hybridization oven. Arrays were washed and stained on an Affymetrix FS450 Fluidics station. The arrays were scanned on an Affymetrix GeneChip Scanner 3000. A summary of the image signal data for every gene interrogated on the array was generated using the Affymetrix Statistical Algorithm MAS 5.0 (GCOS v1.3) algorithm. The cRNA hybridization data from this study are accessible through GEO Series accession number GSE20339 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE20339). Statistical analyses Physiological studies Two-way (accession and treatment) factorial analysis of variance (anova) was used for analysis of all individual soil and plant response variables following tests to assure assumptions of normality and homoscedasticity were met. For these univariate analyses, there were generally data for each variable from six subsample plants for each line and treatment in each of the three replicate blocks; subsample values within each block were averaged before analysis. Soil moisture variables and absolute value of Ysol were ln transformed to meet anova assumptions. All results presented are least-squares (LS) means, back-transformed as needed. Multivariate analysis of water relations responses including leaf RWC, Ytot, Ysol and dC were conducted with multivariate analysis of variance to assess response patterns by Tsu-1 and Kas-1 to the treatments. Statistical analyses of physiological data were conducted with JMP7.0 (SAS Institute Inc., Cary, NC, USA). Comparative genomic hybridizations The Affymetrix ATH1 array was constructed from the full genome sequence of the Col-0 accession of Arabidopsis. Statistical analyses were performed using SAS procedures as called by JMP Genomics or using several procedures implemented in the R statistical package (R Development Core Team 2008). Original CEL files from gDNA hybridizations were imported and processed at the probe level using a Robust Multichip Average (RMA) based background correction, log2 transformation and quantile normalization of raw intensity values (Irizarry et al. 2003). We utilized a custom CDF file (ATH1_AT_TAIR.cdf) constructed from the Arabidopsis Information Resource (TAIR) version 7 of the Arabidopsis genome (available for download at http://brainarray .mbni .med.umich.edu / Brainarray / Database/CustomCDF/.asp) (Dai et al. 2005). This custom CDF file was created using a series of searches to identify unique probes and filter probes with crosshybridization to multiple genomic sites. We confirmed sample labelling by performing a hierarchical clustering of intensity data across all arrays. As expected, replicates of each accession consistently clustered. One-way anova models were subsequently fit for each probe with t-test contrasts between either Tsu-1 or Kas-1 and the Col-0 control line using JMP Genomics 3.2. An empirical Bayes approach was used to shrink the residual variance for each probe based on a prior distribution of the variance estimated from all probes using an invertedgamma distribution. This approach resulted in increased power and sensitivity by improving the stability of the residual variance estimates. We controlled for multiple testing using a positive false discovery rate (pFDR) of 0.05 (Storey 2003). Gene expression analyses We completed analyses of gene expression data using filtered RMA expression values (Irizarry et al. 2003). CEL files were imported into the R environment using the Affy package (Irizarry, Gautier & Cope 2002) and gene expression measures generated using the RMA function (background corrected, log2 transformed, quantile normalized, median-polished summary) with the custom CDF file described above. However, we filtered the processing CDF file of probes that were identified as having significantly different genomic hybridization intensities when compared to the control line Col-0. Filtering was completed with the R package CustomCDF (http://arrayanalysis.mbni. med.umich.edu/MBNIUM.html#CustomCDF), using a 0.05 pFDR criteria, and a requirement of at least three probes per probeset. The removal of probes containing sequence polymorphisms, including single nucleotide polymorphisms (SNPs) and copy number differences, should result in robust gene expression measures that are minimally impacted by sequence divergence among test accessions. We also completed CEL file processing without probe filtering to assess the impact of hybridization polymorphisms on differentially expressed gene lists. Expression measures were subsequently imported into JMP Genomics 3.2 and processed using the anova procedure. In this case, we fit a fixed-effect general linear model including a term for ‘accession’, ‘treatment’ and their interaction using SAS Proc Mixed and an empirical Bayes shrinkage of the residual variance for each probeset. In this case, the shrinkage analysis was completed using the 4 T. E. Juenger et al. © 2010 Blackwell Publishing Ltd, Plant, Cell and Environment method of Ledoit & Wolf (2004) as described by Schafer & Strimmer (2005) with a custom R script. We controlled for multiple testing using a pFDR of 0.05. Quantitative PCR studies In addition to our array studies, we used hybridization based quantitative PCR assays (QPCR) of 31 candidate genes (Table 2) picked from prior studies and a screening of the literature (Bray 2004) to further evaluate gene expression responses. This list contains genes that were found to be consistently regulated by abiotic stress in the microarray studies of Seki et al. (2002) (ABA treatment), Kreps et al. (2002) (mannitol imposed osmotic stress) and Kawaguchi et al. (2004) (progressive soil-water deficit) (Bray 2004) with the addition of several transcription factors of interest (CBF4, DREB2a, DREB1a, DREB1b). Our list is diverse and includes genes involved in metabolism, transport, signal transduction, transcription, hydrophilic/heat-soluble proteins as well as several transcripts of ‘unknown’ function. As described above, rosette leaf material of replicate plants was harvested on RNAlater and extracted as described above for each experimental plant (six biological replicates per accession, per treatment, per block: 6 ¥ 2 ¥ 2 ¥ 3 = 72 plant samples total). We used ProbeLibrary and ABI Taqman assays (genes AT1G01470 and AT2G43570) for our QPCR experiments with ABgene onestep QPCR reagents and the ABI 7900 HT real-time PCR machine. Samples were screened with duplicate technical replicates and subsequent analyses were completed on the replicate averages. Relative mRNA abundance was determined on the basis of the threshold cycle (CT) value for each reaction. We utilized three reference genes (ACTIN, AT3G18780; SAND, AT2G28390; TIP41-LIKE, AT4G34270) as endogenous controls to normalize the quantity of input RNA in reactions. The CT-value for each target gene was subtracted from the geometric mean CT-value for the control genes to obtain DCT-values, which were used in subsequent statistical analyses. DCT-values were analysed using fixed-factor anova with Proc Mixed in SAS (Littel et al. 1996) with block, accession, treatment and accession–treatment interaction effects. We sequenced the probes and primers for each QPCR assay in both Tsu-1 and Kas-1 to evaluate the impact of sequence polymorphisms on expression measures. We found SNPs in assays for three genes (AT1G62570, AT2G40000, AT5G06760), and therefore, the results from these should be taken with caution as polymorphisms in these assays may impact expression estimates.